/**
 * @file FullweightingRes.h
 * @brief The implementation of the fullweighting restrict operator.
 * @author XDDDD
 * @version 
 * @date 2021-06-10
 */
#ifndef __PAISLEYPARK__FULLWRIGHTINGRES_H__
#define __PAISLEYPARK__FULLWRIGHTINGRES_H__

#include "RestrictionOperator.h"

template <size_t DIM>
class FullweightingRes : public RestrictionOperator<DIM> {
public:
	void operator()(MultigridVector<DIM>& _v, const size_t _M) override;
};

template <size_t DIM>
void FullweightingRes<DIM>::operator()(MultigridVector<DIM>& _v, const size_t _M){
	if(_v.get_N() < _M) {
		std::cout << "The present grids do not lies in the coarest grids" << std::endl;
		std::exit(-1);
	}
	else if(_M <= 1) {
		std::cout << "Can not restrict finer." << std::endl;
		std::exit(-1);
	}
	std::vector<size_t> coord_2h;
	std::vector<size_t> coord_h;
	coord_2h.resize(DIM);
	coord_h.resize(DIM);
	int d[DIM] = {0};
	size_t pole = 0;
	bool is_bound = false;
	size_t num_1d = (1 << _M - 1) + 1;
	for(size_t i = 0; i != DIM; i++) {
		coord_2h[i] = 1;
	}
	double coef = 1.0;
	_v.cleanpart(_M - 1);
	//Go through all the nodes in the coareser grid.
	while(pole != DIM) {
		while(pole != DIM) {
				pole = 0;
				coef = 1.0;
				for(size_t i = 0; i != DIM; i++) {
					if(d[i] != 0)
						coef *= 0.5;
					coef *= 0.5;
					coord_h[i] = d[i] +  2*coord_2h[i];
				}
				_v(coord_2h, _M - 1) += coef*_v(coord_h, _M);
				//next
				d[pole] = -2*d[pole] - 1;			//0 -> -1 -> 1 -> -3
				while(d[pole] == -3) {
					d[pole] = 0;
					pole++;
					if(pole == DIM) 
						break;
					d[pole] = -2*d[pole] - 1;
				}
			}
		pole = 0;
		coord_2h[pole]++;
		while(coord_2h[pole] == num_1d - 1) {
			coord_2h[pole] = 1;
			pole++;
			if(pole == DIM)
				break;
			coord_2h[pole]++;
		}
	}
};

#else
//Do nothing.
#endif
